model{
 for(i in 1:n){# n is the number of obs

   xb[i] <- alpha + beta[1]*x[i,1] + beta[2]*x[i,2] + beta[3]*x[i,3] + beta[4]*x[i,4] + beta[5]*x[i,5] + beta[6]*x[i,6] + beta[7]*x[i,7] + beta[8]*x[i,8] + beta[9]*x[i,9] + beta[10]*x[i,10] + beta[11]*x[i,11] + beta[12]*x[i,12] + beta[13]*x[i,13] + beta[14]*x[i,14] + beta[15]*x[i,15]

   lambda[i] <- exp(xb[i])
   p[i] <- rate/(rate + lambda[i])
   y[i] ~ dnegbin(p[i], rate)
 }

  for( j in 1:K){
  beta[j] ~ dnorm(0, 0.1)
  }
  alpha ~ dnorm(0, 0.1)
  rate ~ dunif(0, 1000)
}

